using GeoPhyInv
using Statistics
using Plots; gr();
using ColorSchemes

Medium

medium = Medium(:acou_homo2D, 5); # load a simple homogeneous acoustic medium from the gallery
update!(medium, [:vp, :rho], randn_perc = 5); # add some random noise to the medium
println(medium)
> 2-D Medium with [:z, :x] in [[-1000.0, 1000.0], [-1000.0, 1000.0]]
> number of samples:	[401, 401]
> sampling intervals:	[5.0, 5.0]
> vp:	min	1995.0397603607325	max	3076.4609700370515
> rho:	min	1969.3995953763388	max	3034.211095525826
> Bounds:
5-element Named Vector{Vector{Float64}}
vcat  │
──────┼───────────────────────────
:vp   │           [2250.0, 2750.0]
:rho  │           [2250.0, 2750.0]
:K    │   [1.13906e10, 2.07969e10]
:KI   │ [4.80841e-11, 8.77915e-11]
:rhoI │ [0.000363636, 0.000444444]

AGeom

ageom = AGeom(medium.mgrid, :xwell, SSrcs(2)); # load a simple acquisition using `mgrid` of the medium
println(ageom)
> 2-D acquisition with 2 supersource(s), [1, 1] source(s) and [10, 10] receiver(s)

SrcWav

tgrid = range(0.0, stop = 2.0, length = 2500); # generate a time grid
wav = ricker(15.0, tgrid, tpeak = 0.25); # ricker wavelet
srcwav = SrcWav(tgrid, ageom, [:p]);
update!(srcwav, [:p], wav);

Plotting the vp model

p1=heatmap(medium, :vp, seriescolor=cgrad(:roma))
scatter!(ageom, SSrcs())
scatter!(ageom, Recs())
plot(p1)

Acoustic

Now we have all the required variables to create SeisForwExpt object and prepare the forward modeling. While creating, we switched the snaps_flag on, and instruct recording field at tsnaps. Once the Expt object is created, do the modeling without additional any memory allocations using update!

tsnaps=tgrid[1:div(length(tgrid),20):end] # store 20 snapshots
pa = SeisForwExpt(
    FdtdAcoustic(),
    medium = medium,
    ageom = ageom,
    srcwav = srcwav,
    snaps_field = :p,
    tsnaps = tsnaps,
    tgrid = tgrid,
    rfields = [:p],
    verbose = true,
);
[ Info: frequency bounds for propagating wavefield 1 are: [5.00e-01, 4.10e+01], with peak at: 1.05e+01
[ Info: spatial sampling (5.00e+00) can be as high as 5.49e+00
[ Info: time sampling (8.00e-04) can be as high as 1.29e-03

Update pa to perform modelling; prints time and allocations

@time update!(pa);

FdtdAcoustic() supershot 1/2   1%|▏                     |  ETA: 0:02:18
FdtdAcoustic() supershot 1/2  23%|█████                 |  ETA: 0:00:20
FdtdAcoustic() supershot 1/2  29%|██████▍               |  ETA: 0:00:18
FdtdAcoustic() supershot 1/2  34%|███████▍              |  ETA: 0:00:16
FdtdAcoustic() supershot 1/2  38%|████████▌             |  ETA: 0:00:15
FdtdAcoustic() supershot 1/2  43%|█████████▌            |  ETA: 0:00:14
FdtdAcoustic() supershot 1/2  48%|██████████▋           |  ETA: 0:00:12
FdtdAcoustic() supershot 1/2  53%|███████████▊          |  ETA: 0:00:11
FdtdAcoustic() supershot 1/2  58%|████████████▉         |  ETA: 0:00:09
FdtdAcoustic() supershot 1/2  64%|██████████████        |  ETA: 0:00:08
FdtdAcoustic() supershot 1/2  69%|███████████████▏      |  ETA: 0:00:07
FdtdAcoustic() supershot 1/2  74%|████████████████▍     |  ETA: 0:00:06
FdtdAcoustic() supershot 1/2  79%|█████████████████▌    |  ETA: 0:00:04
FdtdAcoustic() supershot 1/2  85%|██████████████████▋   |  ETA: 0:00:03
FdtdAcoustic() supershot 1/2  90%|███████████████████▊  |  ETA: 0:00:02
FdtdAcoustic() supershot 1/2  95%|█████████████████████ |  ETA: 0:00:01
FdtdAcoustic() supershot 1/2 100%|██████████████████████| Time: 0:00:21

FdtdAcoustic() supershot 2/2   5%|█▏                    |  ETA: 0:00:18
FdtdAcoustic() supershot 2/2  11%|██▍                   |  ETA: 0:00:17
FdtdAcoustic() supershot 2/2  16%|███▍                  |  ETA: 0:00:17
FdtdAcoustic() supershot 2/2  20%|████▌                 |  ETA: 0:00:16
FdtdAcoustic() supershot 2/2  25%|█████▋                |  ETA: 0:00:15
FdtdAcoustic() supershot 2/2  30%|██████▋               |  ETA: 0:00:14
FdtdAcoustic() supershot 2/2  35%|███████▊              |  ETA: 0:00:13
FdtdAcoustic() supershot 2/2  40%|████████▊             |  ETA: 0:00:12
FdtdAcoustic() supershot 2/2  45%|█████████▊            |  ETA: 0:00:11
FdtdAcoustic() supershot 2/2  49%|██████████▉           |  ETA: 0:00:10
FdtdAcoustic() supershot 2/2  55%|████████████          |  ETA: 0:00:09
FdtdAcoustic() supershot 2/2  60%|█████████████▏        |  ETA: 0:00:08
FdtdAcoustic() supershot 2/2  65%|██████████████▍       |  ETA: 0:00:07
FdtdAcoustic() supershot 2/2  71%|███████████████▌      |  ETA: 0:00:06
FdtdAcoustic() supershot 2/2  76%|████████████████▊     |  ETA: 0:00:05
FdtdAcoustic() supershot 2/2  81%|█████████████████▉    |  ETA: 0:00:04
FdtdAcoustic() supershot 2/2  87%|███████████████████   |  ETA: 0:00:03
FdtdAcoustic() supershot 2/2  92%|████████████████████▏ |  ETA: 0:00:02
FdtdAcoustic() supershot 2/2  97%|█████████████████████▍|  ETA: 0:00:01
FdtdAcoustic() supershot 2/2 100%|██████████████████████| Time: 0:00:19
timer output of methods inside fdtd time loop
 ──────────────────────────────────────────────────────────────────────────────
                                       Time                   Allocations
                               ──────────────────────   ───────────────────────
       Tot / % measured:            48.1s / 83.5%            720MiB / 12.0%

 Section               ncalls     time   %tot     avg     alloc   %tot      avg
 ──────────────────────────────────────────────────────────────────────────────
 advance time step      5.00k    37.9s  94.4%  7.58ms   70.6MiB  81.9%  14.4KiB
 record at receivers    5.00k    2.12s  5.27%   423μs   5.74MiB  6.67%  1.18KiB
 add source             5.00k    132ms  0.33%  26.4μs   9.81MiB  11.4%  2.01KiB
 ──────────────────────────────────────────────────────────────────────────────
 48.390481 seconds (12.70 M allocations: 761.451 MiB, 0.61% gc time, 18.07% compilation time)

Extracting snaps

snaps1 = pa[:snaps, 1]; # extracting snaps of the first supersource
snaps2 = pa[:snaps, 2]; # second supersource

Plotting pressure snapshots after acoustic simulation

function plot_snapshots(name)
    clip_perc = 90
    pmax = maximum([maximum(pa[:snaps][it]) for it = 1:length(tsnaps)])
    pmax -= pmax * 0.01 * clip_perc

    anim = @animate for it = 1:length(tsnaps)
        plot(
            [
                (
                    f = pa[:snaps, is][it];
                    heatmap(
                        f,
                        aspect_ratio = 1,
                        xlims = (1, size(f, 2)),
                        ylims = (1, size(f, 1)),
                        yflip = true,
                        seriescolor = cgrad(:seismic),
                        clims = (-pmax, pmax),
                        legend = false,
                    )
                ) for is = 1:2
            ]...,
            title = string("snapshot at: ", round(tsnaps[it], digits = 5), " s"),
        )
    end
    return gif(anim, string(name,".gif"), fps = 1)
end
plot_snapshots("acoustic_snaps")

Elastic

Similarly, we can do the elastic wave-equation modeling

medium = Medium(:elastic_homo2D, 5); # load a simple homogeneous elastic medium from the gallery
update!(medium, [:vp, :vs, :rho], randn_perc = 5); # add some random noise to the medium
println(medium)
ageom = AGeom(medium.mgrid, :xwell, SSrcs(2)); # load a simple acquisition using `mgrid` of the medium
println(ageom)
> 2-D Medium with [:z, :x] in [[-1000.0, 1000.0], [-1000.0, 1000.0]]
> number of samples:	[401, 401]
> sampling intervals:	[5.0, 5.0]
> vp:	min	1963.116052503446	max	3101.0698069092305
> vs:	min	1154.7950463919065	max	1901.4136541003995
> rho:	min	1949.1462081580023	max	3071.836020948651
> Bounds:
8-element Named Vector{Vector{Float64}}
vcat  │
──────┼───────────────────────────
:vp   │           [2250.0, 2750.0]
:vs   │           [1350.0, 1650.0]
:rho  │           [2250.0, 2750.0]
:K    │    [5.92312e9, 1.08144e10]
:KI   │  [9.24695e-11, 1.6883e-10]
:rhoI │ [0.000363636, 0.000444444]
:mu   │     [4.10062e9, 7.48688e9]
:muI  │ [1.33567e-10, 2.43865e-10]

> 2-D acquisition with 2 supersource(s), [1, 1] source(s) and [10, 10] receiver(s)

Add source term on :vz grid

srcwav = SrcWav(tgrid, ageom, [:vz]);
update!(srcwav, [:vz], wav);

Perform modelling

pa = SeisForwExpt(
    FdtdElastic(),
    medium = medium,
    ageom = ageom,
    srcwav = srcwav,
    snaps_field = :vz,
    tsnaps = tsnaps,
    rfields = [:vz],
    tgrid = tgrid,
    verbose = true,
)
@time update!(pa);

snaps = pa[:snaps, 1];
snaps = pa[:snaps, 2];
[ Info: frequency bounds for propagating wavefield 1 are: [5.00e-01, 4.10e+01], with peak at: 1.05e+01
┌ Warning: decrease maximum spatial sampling (5.00e+00) below 2.82e+00
└ @ GeoPhyInv ~/work/GeoPhyInv.jl/GeoPhyInv.jl/src/fdtd/stability.jl:38
[ Info: time sampling (8.00e-04) can be as high as 1.10e-03

FdtdElastic() supershot 1/2   2%|▍                      |  ETA: 0:01:02
FdtdElastic() supershot 1/2   4%|█                      |  ETA: 0:00:46
FdtdElastic() supershot 1/2   7%|█▌                     |  ETA: 0:00:42
FdtdElastic() supershot 1/2   9%|██▏                    |  ETA: 0:00:40
FdtdElastic() supershot 1/2  12%|██▋                    |  ETA: 0:00:39
FdtdElastic() supershot 1/2  14%|███▎                   |  ETA: 0:00:37
FdtdElastic() supershot 1/2  16%|███▊                   |  ETA: 0:00:36
FdtdElastic() supershot 1/2  19%|████▍                  |  ETA: 0:00:35
FdtdElastic() supershot 1/2  21%|████▉                  |  ETA: 0:00:34
FdtdElastic() supershot 1/2  24%|█████▍                 |  ETA: 0:00:33
FdtdElastic() supershot 1/2  26%|█████▉                 |  ETA: 0:00:32
FdtdElastic() supershot 1/2  28%|██████▌                |  ETA: 0:00:31
FdtdElastic() supershot 1/2  30%|███████                |  ETA: 0:00:30
FdtdElastic() supershot 1/2  33%|███████▌               |  ETA: 0:00:30
FdtdElastic() supershot 1/2  35%|████████               |  ETA: 0:00:29
FdtdElastic() supershot 1/2  37%|████████▌              |  ETA: 0:00:28
FdtdElastic() supershot 1/2  39%|█████████              |  ETA: 0:00:27
FdtdElastic() supershot 1/2  41%|█████████▌             |  ETA: 0:00:26
FdtdElastic() supershot 1/2  43%|██████████             |  ETA: 0:00:25
FdtdElastic() supershot 1/2  46%|██████████▌            |  ETA: 0:00:24
FdtdElastic() supershot 1/2  48%|███████████▏           |  ETA: 0:00:23
FdtdElastic() supershot 1/2  51%|███████████▋           |  ETA: 0:00:22
FdtdElastic() supershot 1/2  53%|████████████▏          |  ETA: 0:00:21
FdtdElastic() supershot 1/2  55%|████████████▊          |  ETA: 0:00:20
FdtdElastic() supershot 1/2  58%|█████████████▍         |  ETA: 0:00:18
FdtdElastic() supershot 1/2  61%|█████████████▉         |  ETA: 0:00:17
FdtdElastic() supershot 1/2  63%|██████████████▌        |  ETA: 0:00:16
FdtdElastic() supershot 1/2  66%|███████████████▏       |  ETA: 0:00:15
FdtdElastic() supershot 1/2  68%|███████████████▋       |  ETA: 0:00:14
FdtdElastic() supershot 1/2  71%|████████████████▎      |  ETA: 0:00:13
FdtdElastic() supershot 1/2  73%|████████████████▉      |  ETA: 0:00:12
FdtdElastic() supershot 1/2  76%|█████████████████▍     |  ETA: 0:00:10
FdtdElastic() supershot 1/2  78%|██████████████████     |  ETA: 0:00:09
FdtdElastic() supershot 1/2  81%|██████████████████▋    |  ETA: 0:00:08
FdtdElastic() supershot 1/2  83%|███████████████████▏   |  ETA: 0:00:07
FdtdElastic() supershot 1/2  86%|███████████████████▊   |  ETA: 0:00:06
FdtdElastic() supershot 1/2  88%|████████████████████▍  |  ETA: 0:00:05
FdtdElastic() supershot 1/2  91%|████████████████████▉  |  ETA: 0:00:04
FdtdElastic() supershot 1/2  94%|█████████████████████▌ |  ETA: 0:00:03
FdtdElastic() supershot 1/2  96%|██████████████████████▏|  ETA: 0:00:02
FdtdElastic() supershot 1/2  99%|██████████████████████▊|  ETA: 0:00:01
FdtdElastic() supershot 1/2 100%|███████████████████████| Time: 0:00:42

FdtdElastic() supershot 2/2   3%|▋                      |  ETA: 0:00:38
FdtdElastic() supershot 2/2   5%|█▏                     |  ETA: 0:00:38
FdtdElastic() supershot 2/2   8%|█▊                     |  ETA: 0:00:37
FdtdElastic() supershot 2/2  10%|██▎                    |  ETA: 0:00:37
FdtdElastic() supershot 2/2  12%|██▉                    |  ETA: 0:00:36
FdtdElastic() supershot 2/2  15%|███▍                   |  ETA: 0:00:35
FdtdElastic() supershot 2/2  17%|████                   |  ETA: 0:00:34
FdtdElastic() supershot 2/2  19%|████▌                  |  ETA: 0:00:34
FdtdElastic() supershot 2/2  22%|█████                  |  ETA: 0:00:33
FdtdElastic() supershot 2/2  24%|█████▌                 |  ETA: 0:00:32
FdtdElastic() supershot 2/2  26%|██████▏                |  ETA: 0:00:31
FdtdElastic() supershot 2/2  29%|██████▋                |  ETA: 0:00:30
FdtdElastic() supershot 2/2  31%|███████▏               |  ETA: 0:00:29
FdtdElastic() supershot 2/2  33%|███████▋               |  ETA: 0:00:28
FdtdElastic() supershot 2/2  36%|████████▏              |  ETA: 0:00:28
FdtdElastic() supershot 2/2  38%|████████▋              |  ETA: 0:00:27
FdtdElastic() supershot 2/2  40%|█████████▏             |  ETA: 0:00:26
FdtdElastic() supershot 2/2  42%|█████████▊             |  ETA: 0:00:25
FdtdElastic() supershot 2/2  45%|██████████▎            |  ETA: 0:00:24
FdtdElastic() supershot 2/2  47%|██████████▊            |  ETA: 0:00:23
FdtdElastic() supershot 2/2  49%|███████████▍           |  ETA: 0:00:22
FdtdElastic() supershot 2/2  52%|███████████▉           |  ETA: 0:00:21
FdtdElastic() supershot 2/2  54%|████████████▌          |  ETA: 0:00:20
FdtdElastic() supershot 2/2  57%|█████████████▏         |  ETA: 0:00:18
FdtdElastic() supershot 2/2  60%|█████████████▊         |  ETA: 0:00:17
FdtdElastic() supershot 2/2  62%|██████████████▎        |  ETA: 0:00:16
FdtdElastic() supershot 2/2  65%|██████████████▉        |  ETA: 0:00:15
FdtdElastic() supershot 2/2  67%|███████████████▌       |  ETA: 0:00:14
FdtdElastic() supershot 2/2  70%|████████████████▏      |  ETA: 0:00:13
FdtdElastic() supershot 2/2  73%|████████████████▊      |  ETA: 0:00:12
FdtdElastic() supershot 2/2  75%|█████████████████▎     |  ETA: 0:00:10
FdtdElastic() supershot 2/2  78%|█████████████████▉     |  ETA: 0:00:09
FdtdElastic() supershot 2/2  80%|██████████████████▌    |  ETA: 0:00:08
FdtdElastic() supershot 2/2  83%|███████████████████▏   |  ETA: 0:00:07
FdtdElastic() supershot 2/2  86%|███████████████████▋   |  ETA: 0:00:06
FdtdElastic() supershot 2/2  88%|████████████████████▎  |  ETA: 0:00:05
FdtdElastic() supershot 2/2  91%|████████████████████▉  |  ETA: 0:00:04
FdtdElastic() supershot 2/2  93%|█████████████████████▌ |  ETA: 0:00:03
FdtdElastic() supershot 2/2  96%|██████████████████████ |  ETA: 0:00:02
FdtdElastic() supershot 2/2  99%|██████████████████████▋|  ETA: 0:00:01
FdtdElastic() supershot 2/2 100%|███████████████████████| Time: 0:00:41
timer output of methods inside fdtd time loop
 ──────────────────────────────────────────────────────────────────────────────
                                       Time                   Allocations
                               ──────────────────────   ───────────────────────
       Tot / % measured:            85.8s / 96.8%            233MiB / 47.6%

 Section               ncalls     time   %tot     avg     alloc   %tot      avg
 ──────────────────────────────────────────────────────────────────────────────
 advance time step      5.00k    80.9s  97.4%  16.2ms    109MiB  98.4%  22.3KiB
 record at receivers    5.00k    2.09s  2.52%   419μs    922KiB  0.81%     189B
 add source             5.00k   50.7ms  0.06%  10.1μs    922KiB  0.81%     189B
 ──────────────────────────────────────────────────────────────────────────────
 85.972455 seconds (4.13 M allocations: 251.794 MiB, 0.16% gc time, 3.52% compilation time)

Plotting is vs model

p1=heatmap(medium, :vs, seriescolor=cgrad(:roma))
scatter!(ageom, SSrcs())
scatter!(ageom, Recs())
plot(p1)

Plotting snapshots after elastic simulation

plot_snapshots("elastic_snaps") # P and S waves!